rm(list=ls(all.names=TRUE))
set.seed(221269)


#### Loading data 
load("anes2008.RData")


#### Loading libraries
library(lavaan)
library(semTools)


#### Setting ordered variables ####
orderedvars <- c('respect', 'manners', 'obedience', 'behaved')



################################
#### Configural replication ####
################################

#### Configural model as reported in 
#### Perez & Hetherington 2014, table 2: 
#### Thresholds that are not statistically 
#### different from zero are constrained to zero
#### to release degrees of freedom. 
#### These constraints give three extra 
#### degrees of freedom to the model, 
#### which incidentally improves fit indexes

#### I do not reports these results in the paper

ConfigRepl <- 
  '
## Factor loadings
authorit =~ c(NA,NA)*respect
authorit =~ c(NA,NA)*manners
authorit =~ c(NA,NA)*obedience
authorit =~ c(NA,NA)*behaved


## Thresholds
respect | c(NA,NA)*t1
respect | c(NA,NA)*t2

manners | c(NA,NA)*t1
manners | c(NA,NA)*t2

obedience | c(NA,NA)*t1
obedience | c(00,NA)*t2

behaved | c(00,NA)*t1
behaved | c(NA,00)*t2


## Observed variables (variate): Latent variate scale
respect     ~*~ c(1,1)*respect
manners     ~*~ c(1,1)*manners
obedience   ~*~ c(1,1)*obedience
behaved     ~*~ c(1,1)*behaved


#Latent variables: variances
authorit ~~ c(1,1)*authorit


## Observed variables (latent variate): Intercepts
respect    ~ c(0,0)*1
manners    ~ c(0,0)*1
obedience  ~ c(0,0)*1
behaved    ~ c(0,0)*1


#Latent variable: mean
authorit ~ c(0,0)*1
'

ConfigReplfit <- lavaan::lavaan(
  data=anes2008, 
  model=ConfigRepl, 
  group="race", 
  group.label=c("White", "Black"), 
  ordered=orderedvars, 
  meanstructure=TRUE, parameterization="delta", 
  test="scaled.shifted", se="robust.sem", estimator="ULS",
  mimic="lavaan", model.type="cfa", start="simple", 
  # Weights
  sampling.weights="weight", 
  sampling.weights.normalization="group", 
  # Identification constraints
  std.lv=FALSE, orthogonal=FALSE, auto.fix.single=FALSE,
  int.ov.free=TRUE, int.lv.free=TRUE, auto.fix.first=FALSE,
  auto.th=TRUE, auto.delta=TRUE, auto.var=TRUE, 
  auto.cov.lv.x=TRUE, auto.cov.y=TRUE,
  # Categorical estimation
  zero.cell.warn=TRUE, 
  zero.add=c(0.5,0.5), 
  zero.keep.margins=TRUE, 
  # Computational / optimization options
  check.start=TRUE, check.post=TRUE,
  check.gradient=TRUE, check.vcov=TRUE, 
  check.lv.names=FALSE,
  do.fit=TRUE, optim.method="nlminb",
  implied=TRUE, verbose=TRUE)

summary(ConfigReplfit, fit=T, standardized=T)
lavResiduals(ConfigReplfit)
# fitMeasures(ConfigReplfit)
modificationIndices(ConfigReplfit)


#### Configural model with no constraints
#### on thresholds

Config1.noconstr <- 
  '
## Factor loadings
authorit =~ c(NA,NA)*respect
authorit =~ c(NA,NA)*manners
authorit =~ c(NA,NA)*obedience
authorit =~ c(NA,NA)*behaved


## Thresholds
respect | c(NA,NA)*t1
respect | c(NA,NA)*t2

manners | c(NA,NA)*t1
manners | c(NA,NA)*t2

obedience | c(NA,NA)*t1
obedience | c(NA,NA)*t2

behaved | c(NA,NA)*t1
behaved | c(NA,NA)*t2


## Observed variables (variate): Latent variate scale
respect     ~*~ c(1,1)*respect
manners     ~*~ c(1,1)*manners
obedience   ~*~ c(1,1)*obedience
behaved     ~*~ c(1,1)*behaved


#Latent variables: variances
authorit ~~ c(1,1)*authorit


## Observed variables (latent variate): Intercepts
respect    ~ c(0,0)*1
manners    ~ c(0,0)*1
obedience  ~ c(0,0)*1
behaved    ~ c(0,0)*1


#Latent variable: mean
authorit ~ c(0,0)*1
'

Config1fit.noconstr <- lavaan::lavaan(
  data=anes2008, 
  model=Config1.noconstr, 
  group="race", 
  group.label=c("White", "Black"), 
  ordered=orderedvars, 
  meanstructure=TRUE, parameterization="delta", 
  test="scaled.shifted", se="robust.sem", estimator="ULS",
  mimic="lavaan", model.type="cfa", start="simple", 
  # Weights
  sampling.weights="weight", 
  sampling.weights.normalization="group", 
  # Identification constraints
  int.ov.free=TRUE, int.lv.free=TRUE, auto.fix.first=FALSE,
  std.lv=FALSE, orthogonal=FALSE, auto.fix.single=FALSE,
  auto.th=TRUE, auto.delta=TRUE, auto.var=TRUE, 
  auto.cov.lv.x=TRUE, auto.cov.y=TRUE,
  # Categorical estimation
  zero.cell.warn=TRUE, 
  zero.add=c(0.5,0.5), 
  zero.keep.margins=TRUE, 
  # Computational / optimization options
  check.start=TRUE, check.post=TRUE,
  check.gradient=TRUE, check.vcov=TRUE, 
  check.lv.names=FALSE,
  do.fit=TRUE, optim.method="nlminb",
  implied=TRUE, verbose=TRUE)

summary(Config1fit.noconstr, fit=T, standardized=T)
fitMeasures(Config1fit.noconstr, 
            c('chisq.scaled', 'df.scaled', 'pvalue.scaled', 
              'cfi.scaled', 'rmsea.scaled', 'srmr'))
lavResiduals(Config1fit.noconstr)
modificationIndices(Config1fit.noconstr)


#### Comparing nested models ####
#### This likelihood test
#### compares Perez & Hetherington's 
#### configural model (with impose constraints on thresholds)
#### to the model without such conwtraints. 
#### Results show the models are equivalent.

lavTestLRT(ConfigReplfit, Config1fit.noconstr)





##############################
#### Configural, adjusted ####
##############################

#### This version of the configural model 
#### accounts for the correlation between 
#### the residual covariance RESPECT-BEHAVE
#### Configural model reported in the paper
#### The models below are the ones reported.

Config1.cor.res <- 
  '
## Factor loadings
authorit =~ c(NA,NA)*respect
authorit =~ c(NA,NA)*manners
authorit =~ c(NA,NA)*obedience
authorit =~ c(NA,NA)*behaved


## Thresholds
respect | c(NA,NA)*t1
respect | c(NA,NA)*t2

manners | c(NA,NA)*t1
manners | c(NA,NA)*t2

obedience | c(NA,NA)*t1
obedience | c(NA,NA)*t2

behaved | c(NA,NA)*t1
behaved | c(NA,NA)*t2


## Observed variables (variate): Latent variate scale
respect     ~*~ c(1,1)*respect
manners     ~*~ c(1,1)*manners
obedience   ~*~ c(1,1)*obedience
behaved     ~*~ c(1,1)*behaved


#Latent variables: variances
authorit ~~ c(1,1)*authorit


## Observed variables (latent variate): Intercepts
respect    ~ c(0,0)*1
manners    ~ c(0,0)*1
obedience  ~ c(0,0)*1
behaved    ~ c(0,0)*1


#Latent variable: mean
authorit ~ c(0,0)*1


## Residual correlation, freely est.
respect ~~ c(NA,NA)*behaved
'

Config1fit.cor.res <- lavaan::lavaan(
  data=anes2008, 
  model=Config1.cor.res, 
  group="race", 
  group.label=c("White", "Black"), 
  ordered=orderedvars, 
  meanstructure=TRUE, parameterization="delta", 
  test="scaled.shifted", se="robust.sem", estimator="ULS",
  mimic="lavaan", model.type="cfa", start="simple", 
  # Weights
  sampling.weights="weight", 
  sampling.weights.normalization="group",
  # Identification constraints
  int.ov.free=TRUE, int.lv.free=TRUE, auto.fix.first=FALSE,
  std.lv=FALSE, orthogonal=FALSE, auto.fix.single=FALSE,
  auto.th=TRUE, auto.delta=TRUE, auto.var=TRUE, 
  auto.cov.lv.x=TRUE, auto.cov.y=TRUE,
  # Categorical estimation
  zero.cell.warn=TRUE, 
  zero.add=c(0.5,0.5), 
  zero.keep.margins=TRUE, 
  # Computational / optimization options
  check.start=TRUE, check.post=TRUE,
  check.gradient=TRUE, check.vcov=TRUE, 
  check.lv.names=FALSE,
  do.fit=TRUE, optim.method="nlminb",
  implied=TRUE, verbose=TRUE)

summary(Config1fit.cor.res, fit=T, standardized=T)
lavResiduals(Config1fit.cor.res)
# fitMeasures(Config1fit.cor.res)
modificationIndices(Config1fit.cor.res)



##########################
#### Thresh + Loading ####
##########################

LoadThresh1 <- 
  '
## Factor loadings
authorit =~ c(LoRespWh,LoRespBl)*respect 
authorit =~ c(LoManrWh,LoManrBl)*manners 
authorit =~ c(LoObedWh,LoObedBl)*obedience 
authorit =~ c(LoBehaWh,LoBehaBl)*behaved
LoRespWh == LoRespBl
LoManrWh == LoManrBl
LoObedWh == LoObedBl
LoBehaWh == LoBehaBl


## Thresholds
respect | c(ThRespWh1,ThRespBl1)*t1
respect | c(ThRespWh2,ThRespBl2)*t2
ThRespWh1 == ThRespBl1
ThRespWh2 == ThRespBl2

manners | c(ThManrWh1,ThManrBl1)*t1
manners | c(ThManrWh2,ThManrBl2)*t2
ThManrWh1 == ThManrBl1
ThManrWh2 == ThManrBl2

obedience | c(ThObedWh1,ThObedBl1)*t1
obedience | c(ThObedWh2,ThObedBl2)*t2
ThObedWh1 == ThObedBl1
ThObedWh2 == ThObedBl2

behaved | c(ThBehaWh1,ThBehaBl1)*t1
behaved | c(ThBehaWh2,ThBehaBl2)*t2
ThBehaWh1 == ThBehaBl1
ThBehaWh2 == ThBehaBl2


## Observed variables (variate): Latent variate scale
respect     ~*~ c(1,NA)*respect
manners     ~*~ c(1,NA)*manners
obedience   ~*~ c(1,NA)*obedience
behaved     ~*~ c(1,NA)*behaved


#Latent variables: variances
authorit ~~ c(1,NA)*authorit


## Observed variables (latent variate): Intercepts
respect    ~ c(0,NA)*1
manners    ~ c(0,NA)*1
obedience  ~ c(0,NA)*1
behaved    ~ c(0,NA)*1


#Latent variable: mean
authorit ~ c(0,0)*1


## Residual correlation, freely est
respect ~~ c(ResWh,ResBl)*behaved
# ResWh == ResBl
'

LoadThresh1fit <- lavaan::lavaan(
  data=anes2008, 
  model=LoadThresh1, 
  group="race", 
  group.label=c("White", "Black"), 
  ordered=orderedvars, 
  meanstructure=TRUE, parameterization="delta", 
  test="scaled.shifted", se="robust.sem", estimator="ULS",
  mimic="lavaan", model.type="cfa", start="simple", 
  # Weights
  sampling.weights="weight", 
  sampling.weights.normalization="group", 
  # Identification constraints
  int.ov.free=TRUE, int.lv.free=TRUE, auto.fix.first=FALSE,
  std.lv=FALSE, orthogonal=FALSE, auto.fix.single=FALSE,
  auto.th=TRUE, auto.delta=TRUE, auto.var=TRUE, 
  auto.cov.lv.x=TRUE, auto.cov.y=TRUE,
  # Categorical estimation
  zero.cell.warn=TRUE, 
  zero.add=c(0.5,0.5), 
  zero.keep.margins=TRUE, 
  # Computational / optimization options
  check.start=TRUE, check.post=TRUE,
  check.gradient=TRUE, check.vcov=TRUE, 
  check.lv.names=FALSE,
  do.fit=TRUE, optim.method="nlminb",
  implied=TRUE, verbose=TRUE)

summary(LoadThresh1fit, fit=T, standardized=T)
lavResiduals(LoadThresh1fit)
# fitMeasures(LoadThresh1fit)
lavTestScore(LoadThresh1fit, univariate=TRUE, epc=TRUE)



###################
#### Intercept ####
###################

Intercept1 <- 
'
## Factor loadings
authorit =~ c(LoRespWh,LoRespBl)*respect 
authorit =~ c(LoManrWh,LoManrBl)*manners 
authorit =~ c(LoObedWh,LoObedBl)*obedience 
authorit =~ c(LoBehaWh,LoBehaBl)*behaved
LoRespWh == LoRespBl
LoManrWh == LoManrBl
LoObedWh == LoObedBl
LoBehaWh == LoBehaBl


## Thresholds
respect | c(ThRespWh1,ThRespBl1)*t1
respect | c(ThRespWh2,ThRespBl2)*t2
ThRespWh1 == ThRespBl1
ThRespWh2 == ThRespBl2

manners | c(ThManrWh1,ThManrBl1)*t1
manners | c(ThManrWh2,ThManrBl2)*t2
ThManrWh1 == ThManrBl1
ThManrWh2 == ThManrBl2

obedience | c(ThObedWh1,ThObedBl1)*t1
obedience | c(ThObedWh2,ThObedBl2)*t2
ThObedWh1 == ThObedBl1
ThObedWh2 == ThObedBl2

behaved | c(ThBehaWh1,ThBehaBl1)*t1
behaved | c(ThBehaWh2,ThBehaBl2)*t2
ThBehaWh1 == ThBehaBl1
ThBehaWh2 == ThBehaBl2


## Observed variables (variate): Latent variate scale
respect     ~*~ c(1,NA)*respect
manners     ~*~ c(1,NA)*manners
obedience   ~*~ c(1,NA)*obedience
behaved     ~*~ c(1,NA)*behaved


#Latent variables: variances
authorit ~~ c(1,NA)*authorit


## Observed variables (latent variate): Intercepts
respect    ~ c(IntRespWh,IntRespBl)*1
manners    ~ c(IntManrWh,IntManrBl)*1
obedience  ~ c(IntObedWh,IntObedBl)*1
behaved    ~ c(IntBehaWh,IntBehaBl)*1
IntRespWh == 0 
IntRespWh == IntRespBl
IntManrWh == 0
IntManrWh == IntManrBl
IntObedWh == 0
IntObedWh == IntObedBl
IntBehaWh == 0
IntBehaWh == IntBehaBl


#Latent variable: mean
authorit ~ c(0,NA)*1


## Residual correlation, freely est
respect ~~ c(ResWh,ResBl)*behaved
# ResWh == ResBl
'

Intercept1fit <- lavaan::lavaan(
  data=anes2008, 
  model=Intercept1, 
  group="race", 
  group.label=c("White", "Black"), 
  ordered=orderedvars, 
  meanstructure=TRUE, parameterization="delta", 
  test="scaled.shifted", se="robust.sem", estimator="ULS",
  mimic="lavaan", model.type="cfa", start="simple", 
  # Weights
  sampling.weights="weight", 
  sampling.weights.normalization="group", 
  # Identification constraints
  int.ov.free=TRUE, int.lv.free=TRUE, auto.fix.first=FALSE,
  std.lv=FALSE, orthogonal=FALSE, auto.fix.single=FALSE,
  auto.th=TRUE, auto.delta=TRUE, auto.var=TRUE, 
  auto.cov.lv.x=TRUE, auto.cov.y=TRUE,
  # Categorical estimation
  zero.cell.warn=TRUE, 
  zero.add=c(0.5,0.5), 
  zero.keep.margins=TRUE, 
  # Computational / optimization options
  check.start=TRUE, check.post=TRUE,
  check.gradient=TRUE, check.vcov=TRUE, 
  check.lv.names=FALSE,
  do.fit=TRUE, optim.method="nlminb",
  implied=TRUE, verbose=TRUE)

summary(Intercept1fit, fit=T, standardized=T)
lavResiduals(Intercept1fit)
# fitMeasures(Intercept1fit)

lav_test_int <- lavTestScore(Intercept1fit, univariate=TRUE, epc=TRUE)
restrict <- which(lav_test_int$uni$rhs=="IntRespBl")
(epc_IntRespBl <- lavTestScore(Intercept1fit, release=restrict, 
                              univariate=TRUE, epc=TRUE))
# (epc_interest <- subset(epc_IntRespBl$epc, op=="~1"))



#### Release intercept RESPECT

Intercept2 <- 
  '
## Factor loadings
authorit =~ c(LoRespWh,LoRespBl)*respect 
authorit =~ c(LoManrWh,LoManrBl)*manners 
authorit =~ c(LoObedWh,LoObedBl)*obedience 
authorit =~ c(LoBehaWh,LoBehaBl)*behaved
LoRespWh == LoRespBl
LoManrWh == LoManrBl
LoObedWh == LoObedBl
LoBehaWh == LoBehaBl


## Thresholds
respect | c(ThRespWh1,ThRespBl1)*t1
respect | c(ThRespWh2,ThRespBl2)*t2
ThRespWh1 == ThRespBl1
ThRespWh2 == ThRespBl2

manners | c(ThManrWh1,ThManrBl1)*t1
manners | c(ThManrWh2,ThManrBl2)*t2
ThManrWh1 == ThManrBl1
ThManrWh2 == ThManrBl2

obedience | c(ThObedWh1,ThObedBl1)*t1
obedience | c(ThObedWh2,ThObedBl2)*t2
ThObedWh1 == ThObedBl1
ThObedWh2 == ThObedBl2

behaved | c(ThBehaWh1,ThBehaBl1)*t1
behaved | c(ThBehaWh2,ThBehaBl2)*t2
ThBehaWh1 == ThBehaBl1
ThBehaWh2 == ThBehaBl2


## Observed variables (variate): Latent variate scale
respect     ~*~ c(1,NA)*respect
manners     ~*~ c(1,NA)*manners
obedience   ~*~ c(1,NA)*obedience
behaved     ~*~ c(1,NA)*behaved


#Latent variables: variances
authorit ~~ c(1,NA)*authorit


## Observed variables (latent variate): Intercepts
respect    ~ c(IntRespWh,IntRespBl)*1
manners    ~ c(IntManrWh,IntManrBl)*1
obedience  ~ c(IntObedWh,IntObedBl)*1
behaved    ~ c(IntBehaWh,IntBehaBl)*1
IntRespWh == 0 
# IntRespWh == IntRespBl
IntManrWh == 0
IntManrWh == IntManrBl
IntObedWh == 0
IntObedWh == IntObedBl
IntBehaWh == 0
IntBehaWh == IntBehaBl


#Latent variable: mean
authorit ~ c(0,NA)*1


## Residual correlation, freely est
respect ~~ c(ResWh,ResBl)*behaved
# ResWh == ResBl
'

Intercept2fit <- lavaan::lavaan(
  data=anes2008, 
  model=Intercept2, 
  group="race", 
  group.label=c("White", "Black"), 
  ordered=orderedvars, 
  meanstructure=TRUE, parameterization="delta", 
  test="scaled.shifted", se="robust.sem", estimator="ULS",
  mimic="lavaan", model.type="cfa", start="simple", 
  # Weights
  sampling.weights="weight", 
  sampling.weights.normalization="group", 
  # Identification constraints
  int.ov.free=TRUE, int.lv.free=TRUE, auto.fix.first=FALSE,
  std.lv=FALSE, orthogonal=FALSE, auto.fix.single=FALSE,
  auto.th=TRUE, auto.delta=TRUE, auto.var=TRUE, 
  auto.cov.lv.x=TRUE, auto.cov.y=TRUE,
  # Categorical estimation
  zero.cell.warn=TRUE, 
  zero.add=c(0.5,0.5), 
  zero.keep.margins=TRUE, 
  # Computational / optimization options
  check.start=TRUE, check.post=TRUE,
  check.gradient=TRUE, check.vcov=TRUE, 
  check.lv.names=FALSE,
  do.fit=TRUE, optim.method="nlminb",
  implied=TRUE, verbose=TRUE)

summary(Intercept2fit, fit=T, standardized=T)
lavResiduals(Intercept2fit)
# fitMeasures(Intercept2fit)
lavTestScore(Intercept2fit, univariate=TRUE, epc=TRUE)


#### Comparing nested models: Full v. partial intercept ####
lavTestLRT(Intercept1fit, Intercept2fit)
